source('codes/00_import_libraries.R')


#------------------------------------------------------------------------#
# Functions ----
#------------------------------------------------------------------------#

ols_res = function(y_list, xvar, cntrl, h, yr_range){
  
  # fit an ols model to all returns 
  mod_list = lapply(
    y_list,
    fn_fit_ols,
    xvar    = c(xvar, cntrl),
    df      = dt %>% filter(year %in% yr_range)
  )
  
  # compute newey west se
  mod_robust = mapply(
    function(mod, h){coeftest(mod, vcov. = NeweyWest(mod, lag = h, prewhite=FALSE))},
    mod_list,
    h + 2,
    SIMPLIFY = FALSE
  )
  
  # return the adjusted r2
  is_r2 = lapply(mod_list, function(mod) glance(mod)$adj.r.squared) %>% unlist()
  nobs  = lapply(mod_list, function(mod) glance(mod)$nobs) %>% unlist()
  
  # create huxtale from regression
  out = huxreg(
    mod_robust,
    coefs         = xvar,
    error_format  = "({statistic})",
    stars         = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
    statistics    = c(0),
    borders       = 0,
    outer_borders = 0,
    note          = NULL
  )
  
  # add the r2 to the huxtable
  out = out[-1,-1] %>%
    rbind(is_r2*100)
  
  return(t(out))
}



#------------------------------------------------------------------------#
# Specifications ----
#------------------------------------------------------------------------#

year_grid = data.frame(
  beg = c(1927, 1950, 2000),
  end = c(2019, 2019, 2019)
)
period_list = paste0(year_grid$beg, '-', year_grid$end)


mon_hzn = c(1, 3)
row_names = paste0('$h$=', mon_hzn)


# control lists
vol_list = c(paste0('neg_mkt_l',1:2), paste0('vol_mkt_l',1:2))

vix_list = c(paste0('neg_mkt_l',1:2), paste0('vol_mkt_l',1:2), 
             paste0('vix_l',1:2))



#------------------------------------------------------------------------#
# Data ----
#------------------------------------------------------------------------#

sk = read_excel('input/sp 500 daily.xlsx', sheet = 2) %>% 
  mutate(ym = as.yearmon(Date, '%m/%d/%Y')) %>% 
  filter(year(ym) >= 1928) %>% 
  select(ym, price = Close) %>% 
  mutate(ret = price/lag(price) - 1) %>% 
  drop_na() %>% 
  group_by(ym) %>% 
  summarise(
    mu   = mean(ret),
    n    = n(),
    sig2 = 1/(n-1)*sum((ret - mu)^2),
    vol  = sqrt(sig2),
    sk   = n/((n-1)*(n-2))*sum(((ret - mu)/vol)^3),
    sk_neg = ifelse(sk < 0, sk, 0)
  ) %>% 
  ungroup() %>% 
  mutate(
    sk1 = lead(sk),
    neg_sk1 = ifelse(sk1 < 0, sk1, 0)
  ) %>% 
  select(ym, sk1, neg_sk1)



#------------------------------------------------------------------------#
# Data ----
#------------------------------------------------------------------------#

topic_list = readRDS('input/narrative_dicts_v2.rds')[['base']] %>% 
  names() 
x_list = c(topic_list, 'PLS_All_recursive')
row_names = topic_list %>% str_replace_all('_', ' ') %>% tools::toTitleCase() %>% 
  c('PLS')


dt = readRDS('input/combined_data.rds')[['base']] %>%
  left_join(sk, by = 'ym') %>% 
  mutate(
    across(c(vol_mkt0, vol_mkt1, vol_mkt3), ~ sqrt(250/20 * .x)),
    neg_mkt1 = pmin(mkt1, 0)
  )

for (l in 1:2) {
  dt[paste0('neg_mkt_l',l)] = lag(dt$neg_mkt1, l)
  dt[paste0('vol_mkt_l',l)] = lag(dt$vol_mkt1, l)
  dt[paste0('vix_l',l)] = lag(dt$vix1, l)
}



#------------------------------------------------------------------------#
# Results ----
#------------------------------------------------------------------------#

fn_out = function(yvar, zvar, yr_range, panel) {
  output = lapply(x_list, function(x) {
    ols_res(
      y_list = yvar, 
      xvar = x,
      cntrl = zvar,
      h = 6,
      yr_range = yr_range
    )
  }) %>%
    do.call(what = rbind) %>%
    set_number_format(2) %>% 
    
    insert_column(row_names) %>% 
    insert_row(c('', '$\\beta$', '$t$-stat', '$R^2$')) %>% 
    set_number_format(1, everywhere, 0) %>% 
    set_align(1, everywhere, 'center') %>% 
    set_tb_borders(1, everywhere, 0.8) %>%
    set_bottom_border(final(1), everywhere, 0.8) %>% 
    set_all_padding(0) %>% 
    set_escape_contents(FALSE)
  
  quick_xlsx(output, file = glue('output/tables/Table_9{panel}.xlsx'))
  return(output)
}


fn_out(yvar = 'vol_mkt1', zvar = vol_list, yr_range = 1925:2019, panel = 'A')
fn_out(yvar = 'vix1', zvar = vix_list, yr_range = 1990:2019, panel = 'B')
fn_out(yvar = 'neg_sk1', zvar = NULL, yr_range = 1925:2019, panel = 'C')